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^ , Recently, self-sustained oscillations in complex networks consisting of nonoscil- 

(N ■ 

^ I latory nodes (network oscillators) have attracted great interest in diverse nat- 
^ ' ural and social fields. Due to complexity of network behaviors, little is known 
^ so far about the basic structures and fundamental rules underlying the oscil- 

^ I lations, not to mention the principles of how to control it. In this article we 
■ propose a common design principle for oscillations; predict novel and universal 

Branched Circle (BC) structures of oscillatory networks based on this principle; 
and suggest an operable Complexity Reduction Method to reveal the BC struc- 
tures. These ideas are applied to excitable cell networks (including neural cell 
networks), and genomic regulatory networks. Universal BC structures are iden- 
tified clearly in these two considerably different systems. These BC structures 
reveal for the first time both oscillation sources and wave propagation pathways 
of complex networks, and guide us to control the oscillations with surprisingly 
high efficiency. 
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Introduction 



Self-sustained oscillations in complex networks consisting of non-oscillatory nodes are very 
popular phenomena in natural and social systems. These oscillations are extremely important 
in controlling various basic rhythms in wide fields, such as oscillatory neural networks (1-7), 
sinoatrial node rhythms in cardiac systems (8-10), oscillatory cycles in genomic regulations 
(11-17) and so on. Though the topic of self-sustained oscillations has been investigated 
for some decades, many fundamental questions remain in puzzle. For instance,we do not 
know whether there are some common principles underlying the oscillatory behaviors of 
complex networks in diverse fields; whether there are some common structures hidden in 
the complicated interaction schemes that determine the dynamics of oscillations; and if 
these common structures do exist, how one can find them. Specially, an given oscillatory 
networks as Figs. lA, 15, one can hardly say anything about where the oscillation sources 
are; how oscillatory waves propagate from the sources to the whole networks and how one can 
efficiently control the oscillations based on these understandings. None of these questions of 
crucial importance has been answered if networks are sufficiently complicated. 

In this article we have made the following essential advances, (i) We start from a simple 
while solid basis design principle: any individually nonoscillatory node can oscillate if and 
only if it is driven by one or few interactions with advanced phases, (ii) Based on this design 
principle, all periodically oscillatory one-dimensional (ID) networks must have universal 
unidirectional coupling structures of branched circles (BCs, Fig. 2) in the form of circles with 
radiating branches, (iii) Oscillatory high-dimensional complex networks can be reduced to 
ID EC networks by applying the method of dominant phase-advanced driving pathes. 

We further apply the above ideas to both oscillatory excitable cell networks (ECNs, neural 
cell networks included) (5-7,18-20) and genomic regulatory networks (GRNs) (21-24). It's 
the first time the oscillations of these two considerably different systems are studied with 
a unified approach. We have reached the following results, (i) We apply the same design 
principle and dimension reduction approach to two kinds of systems and obtain the same 
universal EC structures; (ii) From the EC patterns of both systems we can clearly identify 
the oscillation sources (unidirectional regulatory circles) and reveal the phase propagation 
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pathways (unidirectional tree branches); (iii) With these BC patterns we can classify one or 
few most important nodes, by controlling which we can control the oscillations of the whole 
networks with surprisingly high efficiency. 

Design principles and universal structures 

Figure lA considers an ECN example in which multiple cells with a large number of random 
couplings show a rather complicated interaction pattern. Although each individual node is 
nonoscillatory, with certain initial conditions we observe periodic oscillations, one of which 
is shown in Fig. IB. We further study how sensitive the oscillation is to the control. After 
exhaustive tests we find that the oscillation can be terminated (i.e., turned to the homoge- 
neous rest state ui = Vi = 0,i = 1,2, .., N) by removing a single red node (node 78) shown in 
Fig. 1 C, while the oscillation persists safely if we remove as many as 70 empty square nodes 
in Fig. IC. Here removing a node we mean to discard all interactions from and toward the 
given node. The results of Fig. IC are highly surprising. They clearly show that though 
the topological interaction structure of Fig. lA is random and homogeneous, the dynamic 
organization supporting self-sustained oscillation is strongly heterogeneous, and the roles 
played by different nodes in this organization are considerably different. The central task 
of the present work is to reveal these self-organized patterns under the conditions of full 
knowledge of coupling structure and output data and then to achieve an effective control of 
the oscillatory networks based on this understanding. 

Considering a network with N nodes, dynamic variables arc associated to each node, and 
these variables obey well defined coupled ordinary differential equations (ODEs). Each node 
is nonoscillatory individually while the entire complex networks are periodically oscillatory. 
Regardless of different dynamics and coupling forms for different systems, we propose a 
common design principle for such oscillatory networks. 

Design principle: each nonoscillatory node can oscillate if and only if it is driven by one 
or few oscillatory interactions with advanced phases. 

The definitions of "advanced phase" for different systems will be explained later. Let us 
first consider the simplest ID oscillatory network with each node phase-advancedly driven 
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by one node only (networks with N nodes and N unidirectional interactions). Suppose an 
arbitrary node ii is phasc-advancedly driven by a node 12 via coupling, which is phase- 
advancedly driven by node is in turn, and this successive unidirectional driving chain goes 
as ii Z2 <— is ijfc . Since N is finite we must come to a node iq, q < N, 

which is driven by one of the previous nodes Zi, ^2, iq-i, say ip {p < q). Then a successive 
regulatory loop ip <— ip+i ig <— ip is formed, serving as the oscillation source of all 

other nodes. Therefore, all ID oscillatory networks must have the structure of circles with 
radiating branches, schematically shown in Fig. 2. 

The pattern of Fig. 2 gives a picture of a branched circle (BC), and it is thus called 
as BC structure. This structure is universal for self-sustained periodic oscillations in ID 
networks consisting of nonoscillatory nodes. Since no nonoscillatory node oscillates without 
phase-advanced driving from other nodes, two key rules must be obeyed by any BC structure: 

(i) There must be few (at least one) successively phase-advanced driving circles. 

(ii) Each node not in the circles must be in a tree branch rooted at a node in a circle. 
The BC circle is obviously the oscillation source without which the network can never 

oscillate, and the tree branches show various pathways of phase propagations starting from 
different circle nodes. All nodes in Fig. 2 can be classified according to their locations in the 
BC pattern. We expect that the circle nodes controlling large branches may be of the most 
importance for the oscillation. In comparison all nodes near the branch ends with few or 
even without downstream nodes have the lowest infiuences on the oscillation of the network. 
We will show later that these expectations are well confirmed by numerical results with large 
probability. 

The simple and instructive scheme of Fig. 2 is deduced in ID phase-advanced driving 
networks. However, interaction structures of complex networks in general (e.g.. Fig. lA 
which are high-dimensional and random) are much more complex than Fig. 2. Therefore, we 
propose an operable and physically meaningful method to reduce original random patterns 
(as Fig. lA) to the simple and instructive BC pattern of Fig. 2. The method consists of the 
following Complexity Reduction (CR) steps: 

(a) Find phase- advanced driving interactions for each node. 

(b) Find the single dominant interaction among these phase-advanced driving interac- 
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tions. 



(c) Use all these dominant interactions to unidirectionally link the network nodes, and 
draw the dominant phase-advanced path pattern, which turns out to be the ID BC pattern 
of Fig. 2. 

All the above steps are generally applicable in diverse fields for self-sustained oscilla- 
tions of complex networks of individually nonoscillatory nodes. The particular meanings 

of "advanced phase" and "dominant phase-advanced driving" should be properly defined, 
according to realistic physical, chemical and biological interaction mechanisms in each indi- 
vidual system. 

Coupled excitable cell networks 

We now consider complex excitable cell networks (ECNs) of the Bar Model(18) 



in which Uij means the variable u of the jth node linked to node i. Note, two major 
features of neural networks arc precisely the excitability of cell dynamics and the complexity 
of interaction network(20, 25-27). Without couplings all cells of ECNs are not oscillatory 
individually for certain given a, b, they evolve asymptotically to the rest state u — v — 
and will stay there forever unless some external force drives them from this state. Therefore, 
all the analyses in the former section are applicable to this type of systems. Whenever 
a cell is kicked from the rest state by a small stimulus, the cell can oscillate by its own 
internal dynamics (so called excitable dynamics). Therefore, for a given node that enters 
the region of the rest state {u < Utu) at time tg and departs from this region at time tg, we 
define "phase-advanced drivings" by the interactions from those neighbors which leave from 




dt 



f{ui) - V, 



i = l,2,...,Ar 



(1) 
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the rest state earlier than the given node (i.e., in the period (ts,te)) which thereby provide 
favorable interactions in kicking the given cell from the rest state. Among these phase- 
advanced interactions the dominant phase- advanced driving is defined by the interaction 
from the node first leaving from the rest state in the period ( tg, te)- It is no doubt that the 
dominant driving must give the most important contribution to excite the given node, i.e., 
to drive the given node to oscillate. 

For simplicity, we assume symmetric couplings, and adopt random complex networks 
with identical coupling degree u (i.e. each cell couples to equal number u of other cells). 
We also take identical parameters for all cells. One advantage of this simplest homogeneous 
assumption is to make sure that all self-sustained oscillatory behaviors here are not due to 
any heterogeneity in topological structure, but due to the self-organized heterogeneity of 
dynamical mutual excitations. 

Figure lA shows an ECN with random mutual couplings with identical degree u — 4. 
With the structure of Fig. lA we simulate the system by taking different sets of random 
initial conditions. In most of cases the system evolves asymptotically to the homogeneous 
rest state. However, about 8% of tests provide periodic oscillations, and the state in Fig. IB 
is one of them. 

In Fig. 3^ we plot the BC pattern draw from Figs. lA and IB hy applying the CR 
method. We find a BC pattern being the type of Fig. 2. In Fig. 3 A the BC circle plays 
the role of oscillation source, and phase waves propagate down all the tree branches rooted 
at various cells in the circle. The BC structure of Fig. 3A clearly shows distinctive levels 
of significances of different cells for oscillation that cannot be observed in Fig. 1^1. In Fig. 
lA all cells stand equivalently in the homogeneous and randomly coupled network, and no 
cell takes any priority over others from the topological structure. The situation in Fig. 3A 
is different. Cells in the circle and cells in various turning points of large branches (which 
control large numbers of downstream cells) are likely to be important to the oscillation. In 
particular, node 78 is likely to be the most important cell because it locates in the circle 
on one hand and controls large branches with a huge number of downstream nodes on the 
other hand. It is interesting to observe that node 78 is the very single red node in Fig. 1 C 
that we found important for controlling the oscillation but did not know the reason then. 
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In Fig. SB we show how the oscillation collapses quickly to the homogeneous rest state after 
removing only a single red node 78. On the other hand, cells near the branch ends may be 
much less significant for the oscillation. We remove simultaneously a large number of cells 
(70 nodes) in branch tails (See Fig. 3(7) which are exactly the empty square nodes in Fig. 
IC, the network not only continues its periodic oscillation (Fig. 3D), but also keeps the BC 
structure almost unchanged for the cells remained (Fig. 3C}. Therefore, the question for Fig. 
IC, why so many empty square nodes together are much less important than a single red 
cell, is clearly answered in Fig. 3A: because all these empty square nodes are far from the 
centered oscillation generator, and have little influence on the self-sustained oscillation; on 
the other hand the single red cell 78 controls the oscillation source and a huge number of 
downstream nodes, and it is of crucial importance for the given oscillation. Removing one 
or few other cells in the circle may not stop the oscillation but can considerably change the 
BC structure of the oscillation source. 

Frequency is an important quantity describing the properties of oscillatory networks. We 
further study the influence of BC circles on frequencies of networks. We computed Eqs. (1) 
by taking different random couplings and random initial conditions for N = 100, ly = 3, and 
found some oscillatory realizations (all are periodic), we then measured the frequency u of 
each oscillatory network, and plot a; vs n in Fig. 3E with n being the size of the corresponding 
circle. The red square represents the average frequency < a; > and the solid line denotes the 
linear fltting of the tendency. In Fig. 3E monotonous and nearly linear decrease of < a; > 
with n is clearly demonstrated. These observations convincingly support the conclusion that 
BC circles play the role of oscillation sources in complex networks. 

Figure. 4:A is another ECN network with — 100 and u = 3. For certain initial condition, 
we observe periodic oscillation of Fig. 45. In Fig. AA we show that this oscillation can be 
terminated by removing a pair of nodes (red nodes 12 and 21) in contrast with Fig. IC where 
oscillation can be suppressed by removing only a single node. Similar to Fig. 1 C oscillation 
persists when we simultaneously remove all 60 empty square nodes shown in Fig. 4A. The 
mystery difference between Figs. 1 C and AA can be again well explained by the corresponding 
BC patterns. In Fig. 4C we show the BC pattern corresponding to the oscillation in Figs. 
A A and 45. The essential and interesting difference between Fig. 3 A and Fig. 4Cis that Fig. 
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4C contains two BC clusters instead of the single one in Fig. 3 A. The two red nodes shown 
in the two BC circles of Fig. 4C are exactly identical to those in Fig. 4A with the same color. 
Comparing the BC patterns of Fig. 4C with Fig. 3A, we understand the difference of Fig. 
IC and Fig. A A immediately. Since there are two oscillation source circles in Fig. 4C, we 
have to destroy both circle structures for terminating oscillation by removing two key nodes 
in Fig. 4:A simultaneously, each from a circle of Fig. 4C, and this is sharply different from 
the single circle structure of Fig. 3A. In Figs. AD and AE we get the BC patterns with one of 
the two red nodes (node 12 and node 21) removed, respectively. It's interesting to see that 
in both cases one BC circle of Fig. 4C is destroyed while the other remained circle serves 
as the unique oscillation source, and all nodes of the destroyed BC cluster are engrafted to 
the survival circle for continuing their periodic oscillations. In Fig. AF we remove 60 side 
cells (the empty nodes) simultaneously from Fig. AC, and find again that periodic oscillation 
persists, and the corresponding BC pattern of the unremoved nodes are not affected. We 
have tested a number of different ECNs with different initial conditions, different interaction 
degree and structures, different system sizes, or different ECN models including FHN neural 
cell networks, and obtained very rich behaviors of BC patterns. The general conclusions 
of universal BC patterns are verified in all cases, and some of these results are shown in 
Supporting Information Part 1. 



Complex genomic regulatory circuits 

We now consider another model of self-sustained oscillations of genomic regulatory networks 
(CRN), the dynamics is described by the following coupled ODEs(21-24). 



dx. 
Itt 



Ai{x) Positive regulation , 

Ri{x) Negative regulation , 

Ai{x)Ri{x) Joint regulation , 



X = {xi,X2, ■ ■ .,Xn) , 



Ai{x) = 



act 



hi 



acp + K 



^, Ri{x) ^ {I - ^li)- 



K 



rep'l' + K'- 



(2) 



8 



N N 

acti = XI ' ^^Pi = Yl l^f^i ' i = 1' • • • ' ^) ' 

where xi represents the concentration of protein corresponding to node i, and acti ijepi) 
represents the summation of activatory (repressive) transcriptional factors. These ODEs can 
be derived from a full set of equations of both mRNA and potein concentrations via adiabatic 
approximation when the time scales of transcription and translation are separable(ll). For 
simplicity we consider homogeneous parameter distributions again for all nodes {jii — 
'ji = J, hi = h, Ki = K, ai = /3i = 1, i = 1,2, N). It's emphasized that the GRN 
dynamics Eqs. (2) is apparently different from ECN Eqs. (1) for their intrinsic node dynamics 
(one-variable passive dynamics for Eqs. (2) against two-variable excitable dynamics for Eqs. 

(1) ), coupling dynamics (highly nonlinear positive or negative regulatory interactions for Eqs. 

(2) while linear and diffusive coupling for Eqs. (1)), and coupling structure (unidirectional 
couplings for Eqs. (2) against symmetric couplings for Eqs. (1)). It is a surprise for us 
to observe essentially the same dynamic BC structures in both GRN and ECN systems as 
shown in the following. 

Each node in Eqs. (2) has passive dynamics. Without coupling, variable Xi must evolve 
to a fixed value, and any periodic oscillation of Xi must be driven by one or few periodic 
interactions from other nodes. Let us approximately simplify an arbitrary one-variable pas- 
sive dynamics with a periodical driving x = /j, — jx + f{x,t) by linearizing the oscillatory 
elements around a stable stationary solution of the autonomous dynamics 



Ax ^ -X Ax + Acos{ujt) , (3) 
which has the asymptotic periodic motion 



VA2 + a; 

sin(t) = .^^ ^ , = J" , (4) 



+ a;2 - ^prp 



leading to the phase- advanced driving condition of 
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O<0<f- (5) 

We represent the phase of node i by 0^ and the phase of interaction from node j to node i 
by (pj^i, which is identified by 

/■ 

positive regulation , 
= I (6) 
(t)j + TT negative regulation . 
The condition of phase-advanced driving reads 

0<4>j^i-4>i<l . (7) 

With interaction structures and full periodic oscillation data of complex GRN known we 
can use Eqs. (6) and (7) to identify all phase-advanced interactions among all interactions 
(CR method (a)). To operate the CR method (b) we simply define the dominant driving 
by the interaction among all phase-advanced interactions which has the largest oscillation 
amplitude. The detailed discussion on the phase and amplitude computation is given in 
Supporting Information Part 2. After choosing the dominant interactions for various nodes, 
we can link all nodes by the dominant interactions and construct BC patterns from the 
original GRNs. 

It is known that a necessary condition for a genomic regulatory loop to be oscillatory is 
that the genes in the loop must interact successively in a manner of negative feedback(13,21- 
24), i.e., the number of negative couplings should be odd(13,22). We call a successive uni- 
directional interaction loop with an odd number of negative interactions as an oscillatory 
negative feedback loop(NFL)(28) in the following. Note, existence of one or multiple oscilla- 
tory NFLs is the necessary but not sufficient condition of oscillatory networks. In Supporting 
Information Part 2 (Fig. S5), we show these oscillatory NFLs for q ^ 5. Nevertheless, given a 
complex network (e.g.. Fig. 5D or 5E), there may exist a huge number of possible oscillatory 
NFLs. There is so far no report of a method to find out which NFLs are in function to 
produce a given oscillatory pattern. And this is right our following task. 
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In order to examine the validity of the general prediction on universal oscillatory BC 
patterns wc have computed huge number of GRNs by varying dynamic parameters (homo- 
geneous and heterogeneous), changing autoregulatory dynamics (with or without autoreg- 
ulation, with negative or positive autoregulation) ; changing the number of nodes and the 
structure of couplings and also by varying the initial variable preparations. Specifically, we 
have made 328400 tests (about 1400 oscillatory reahzations) by exhaustively computing dif- 
ferent coupling structures of 3-node circuits; iC random tests for 4-node and 5-node circuits 
(totally 2 X 10^ tests, 708 oscillations); 5 x 10^ random tests for 10-node and 20-node circuits 
(totally 1 X 10^ tests, about 1083 oscillations), we found that none of the tests violates the 
predictions. In particular, we found that in all oscillatory cases we can succeed in construct- 
ing BC patterns all of which have circles being one of the oscillatory NFLs in Fig. S5. In Figs. 
5A-E we show some periodically oscillatory examples of regulatory networks. Although the 
complicated interaction structures don't disclose any clue of the mechanism supporting the 
oscillations, by applying CR method we succeed in reducing the original complex networks 
to the corresponding BC patterns greatly simplified in Figs. 5F-J, respectively, which fully 
confirm the prediction of Fig. 2. Each BC pattern in Fig. 5 has a source circle being one of 
the oscillatory NFLs in Fig. S5, and all other nodes are in various tree branches rooted at 
one of circle nodes, showing wave propagation pathways. 

From Figs. 5F-J we expect that the nodes in the circles or near the circles may be 
important for the given oscillatory states while nodes near the ends of various branches may 
be less significant. We study each gene's influence on oscillations by removing it. In Figs. 5F- 
J any node whose individual removal results in the termination of the oscillations are filled 
with red color, and empty otherwise. We find that by removing a gene on a circle we have 
very large probability to terminate the oscillation. However, when we remove a node located 
at the end of a branch pathway oscillations have much larger probability to persist. For 
statistics we have made a detailed investigation for 10-node oscillatory GRNs, and found that 
if we remove an arbitrary single node on BC circles the probability to terminate oscillations 
is about 84% while this probability is down to about 24% if an arbitrary single node on 
branches is removed. 

For identifying the system response to control, we study the dynamic behavior oi N — 20 
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(Fig. 57) in more detail. In Fig. 6 A we show oscillation of < x{t) >= Xi{t) damping 

to a fixed value after a key circle node removed at t = 1000. On the other hand, the 
oscillation persists (Fig. 6B) and the BC structure is only slightly modified (Fig. 6(7) as a 
node at a branch end is removed. Most interestingly, whenever self-sustained oscillations 
are maintained after removing some nodes, the BC circles have a strong tendency to be 
unchanged (Fig. 6C^ or slightly modified by refinding some interaction bridges to repair the 
circles (Fig. 6D). All these observations verify the significance of the universal BC structures 
for oscillations of complex networks. 

In Eqs. (2) we assume "AND" role between activators and repressors for multiple-factor 
regulations(21). Some regulatory circuits may obey "OR" role(ll,23,29). Though the cou- 
pling dynamics of "OR" rule looks considerably different from Eqs. (2), all analyses for Eqs. 
(2) can be identically applied to the "OR" cases. This aspect is discussed in Supporting 
Information Part 3. 

Discussion 

In conclusion we study the problem of self-sustained periodic oscillations in complex net- 
works consisting of nonoscillatory nodes. We propose a general design principle of oscillatory 
networks, based on which we reveal that phase-advanced driving BC patterns (Fig. 2) are the 
universal structures of simplest ID oscillatory networks. And complicated high-dimensional 
networks can be reduced to these ID BC patterns by applying the method of dominant 
phase-advanced interactions. From the BC patterns we can easily identify oscillation sources 
and phase propagation pathways of oscillatory complex networks. All these messages are 
deeply hidden in the original complex coupling structures and random phase distributions. 
These BC structures are extremely important for understanding and efficiently controlling 
self-sustained oscillations of complex systems. We successfully used these ideas and methods 
to analyze models of excitable cell networks and genomic regulatory circuits. These ideas, 
methods and universality of structures are expected to be applicable to self-sustained os- 
cillations of complex networks in broad range of fields. In recent decades, the concept and 
functions of central pattern generators (CPGs) have attracted great attention in the field of 
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neural networks(30-33). In this article we show for the first time how to uncover CPG-hke 
patterns in periodically oscillatory complex networks from complicated interaction structures 
and seemingly mess phase distribution data. 

In the present article we consider only cases of periodic oscillations where all BC patterns 
are stationary. If oscillations are quasiperiodic or even chaotic, BC patterns may vary during 
the evolutions, and this opens a new field for the further study. Moreover, throughout this 
article we study how to reveal BC patterns with full knowledge of the interaction structures 
and oscillation data. These conditions are not fulfilled in many experiments. Thus, it is 
another crucial task to extend the investigations to the cases with partial data available. 

This work was supported by the National Natural Science Foundation of China un- 
der Grant Nos. 10675020 and the National Basic Research Program of China (973 Pro- 
gram) (2007CB814800) . 
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Fig. 1. A complex excitable cell network of Eqs. (1) consisting of N — 100 nodes, a — 0.84, 
b = 0.07, s = 0.04. All these parameters will be used in Figs. 1, 3, 4. All couplings with 
strength D„ = 0.2 are randomly chosen between N nodes each having coupling degree u — 4. 
(A) An example of network interaction structure. (B) Periodic orbit of average variables 
< u{t) >— j^'^fLiU'iit) and < v{t) >— jf'^f=-iVi{t) of the network (A) for a set of initial 
conditions randomly chosen. (C) Oscillation of (B) can be terminated by removing a single 
red node. By removing a node we mean to discard all interactions from and towards the 
given node. However, oscillation persists if we simultaneously remove other 70 empty square 
nodes. 

Fig. 2. Schema of universal structure of self-sustained oscillatory ID networks with all nodes 
nonoscillatory individually. The arrowed lines indicate unidirectional phase- advanced inter- 
actions. Characteristic features are: one or few unidirectionally interacting circles serving as 
the oscillation sources together with unidirectionally interacting branches radiated from the 
circle showing wave propagation pathways (called branched circles, EC). 
Fig. 3. EC patterns and their oscillation dynamics of ECNs. (A) EC pattern of Eqs. (1) 
reduced from Figs. 1(A) and 1(B) by applying CR method. From this pattern we are able 
to identify oscillation source (the unidirectional circle) and wave propagation pathways (the 
tree-like branches from various nodes of the circle). Self-sustained oscillation of Fig. 1(B) 
can be effectively suppressed by removing only a single red node 78, which is the same as 
that in Fig. 1(C). (B) Oscillation suppression by removing node 78 at t — to — 10. (C) EC 
pattern drawn after removing simultaneously the 70 empty square nodes shown in Fig. 1(C). 
Oscillation is kept and the BC pattern of the remaining nodes is not affected. (D) Oscillation 
evolution with 70 empty square nodes of (C) removed at t = t£) = 10, simultaneously. (E) 
Frequency a; of a periodic complex network of Eqs. (1) plotted vs the size n of the EC 
circle. Large red square represents the average frequency < cu >, and solid line is the linear 
fitting of the data. N = 100, u = 3, = 1.0. Coupling structures and initial conditions 
are randomly chosen. A strong correlation between < u > and n is demonstrated by the 
monotonously decreasing curve. 

Fig. 4. Manipulations of oscillatory complex networks. (A) Another ECN network with 
N — 100, u — 3, Du — 1.0. (B) Oscillatory orbits of network (A) with a certain set of initial 
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condition. Oscillation persists after all the 60 empty square nodes shown in (A) are removed 
simultaneously. Apparently different from Fig. 1(C), now removing any single node can no 
longer terminate the oscillation. Oscillation can be suppressed by, at least, removing the 
pair of red nodes in (A). The difference between Figs. 1(C) and 4(A) is again an interesting 
mystery. (C) BC pattern of the oscillation in (B) considered. Now we find two separated 
BC clusters from two circles. Removing the pair of red nodes 12 and 21 can break the two 
source circle and terminate the oscillation. (D) BC pattern with a single node 12 removed. 
Now one BC circle with node 12 is destroyed and the remained circle of (C) plays the role of 
the unique oscillation source. (E) BC pattern with node 21 removed. (F) BC pattern after 
removing 60 side nodes from (C) (i.e., the 60 empty square nodes in (A)). 
Fig. 5. Complex oscillatory GRNs and the corresponding BC patterns. The sohd green 
(red dashed) line denotes positive (negative) regulation. The lines in the following figures 
on GRNs have the same meanings throughout the article. (A)-(E) Some examples of 4- 
node (/ = 9 interactions), 6-node (7 = 20), 10-node (/ = 55), 20-node (/ = 105), 40-node 
(/ = 224) GRNs , respectively. //^ = = 0, 7^ = 7 = 0.1, hi ^ h ^ 2, Ki ^ K ^ 0.3, 
cti = Pi = 1, i = 1,2, ... ,N. With the given interaction structures all these circuits show 
self-sustained periodic oscillations with arbitrary initial conditions. (F)-(J) BC patterns 
constructed from the periodic oscillation of networks of (A)-(E), respectively, by applying 
CR method. In all cases we find that each BC circle is one of the oscillatory NFLs in Fig. 
S5. A node whose individual removal terminates the oscillation is filled with red, otherwise 
it is empty. 

Fig. 6. Detailed numerical results of Eqs. (2) for the case of Figs. 5(D) and 5(1) (A^ = 20). 

(A) Time evolution of the average concentration < x >, which periodically oscillates for 
t ^to and collapses to a stationary solution after we remove a circle node 15 at to — 1000. 

(B) (C) Time evolution and BC pattern , respectively, with node 6 removed at to — 1000. 
The BC circle of Fig. 5(1) is not changed by removing this branch gene. (D) BC pattern 
with node 18 removed. In (D) the BC circle must be changed by removing a circle gene, 
however, the system repairs the broken circle by finding some interaction bridges with most 
of original circle nodes kept on the new circle. 
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Supplementary Information 



Part 1 Oscillatory excitable cell networks with random 
symmetric interactions 

We have investigated a large number of different ECN. First, we liave tested ECNs of Eqs.(l) 
for different random coupling structures, different random initial conditions and different sys- 
tem sizes. Prom these tests we found many oscillatory configurations. For each oscillatory 
network we drew the dominant phase-advanced interaction pattern, and found that all pat- 
terns show BC structures of type Fig. 2. 

In Fig. Sl^ we show an interaction structure of an ECN with N — 200, v — 3 and 
Du — 1.0, which is periodically oscillatory for certain initial conditions. Though the random 
interaction network looks even more complicated than Fig. lA, the reduced BC patterns 
are still simple and instructive. In Figs. SIB-D we show three different BC patterns for 
three given initial conditions. These oscillations can be terminated by removing only one 
or two key red circle nodes (e. g. node 30 for Fig. SIB, nodes 6 and 38 for Fig. SIC; and 
nodes 126 and 178 for Fig. SID). A particular case is Fig. SID where the BC pattern 
has a single long An circle consisting of 21 cells (by 2kn circle we mean that phase angle 
variation in the circle is 2kn with wave number k). Now wc can't suppress oscillations by 
removing a single cell. Instead, a pair of cells with phase angle distance about 27r should be 
removed simultaneously. In all the three cases of Figs. SlB-D oscillations can persist and BC 
circles cannot be slightly changed after 70% nodes (all at the end of branches) are removed 
simultaneously. One of these results is given in Fig. Sl£^. It is noticed that in Figs. SIB-D the 
nodes immediately upstream to each red node, provide the only phase-advanced interactions 
for the corresponding red nodes. Therefore, removing the nearest upstream node of any red 
node is equivalent to remove this downstream red node. 

In Fig.S2 we show interesting BC pattern manipulations when oscillations are not termi- 
nated by removing some key nodes. In Fig. S2A{S2B) we remove a red cell 6 (38) from Fig. 
SI C, and find that the up-left (low-right ) BC circle is destroyed while the other remaining 
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circle serves as the only source circle in the new BC pattern. In Fig. S2C {S2D) we remove 
a single red node 126 (red node 178) from Fig. SID. Now we find that the BC circle shrinks 
from the 47r circle to a 27r loop with the key node 178 (key node 126) serving as the key 
circle node in the new BC structure. Pictures of Figs. 4 and S2 show various structure 
reorganizations by controlling few key nodes recognized from the BC patterns. On the other 
hand, oscillations are insensitive to the control of various branch nodes. 

In Fig. S3 we plot various BC patterns for networks with even larger sizes. In Figs. S3A'C 
we take A'" — 400 for different interaction structures and different initial conditions. We find 
BC patterns with a 27r circle (A), An circle (5), and two 27r circles (C), respectively. The 
numbers inside the circles indicate the indexes of the circle nodes while those outside the 
circles associated with arrows denote the size of branches rooted at the given circle nodes. 
In Figs. S3 D-F we do the same as Figs. S3A-Chj taking = 2500. Now the network size 
is many times larger than Fig. lA. However, the prediction of the universal BC structure of 
Fig. 2 is still verified perfectly. An interesting phenomenon is that we find triple BC clusters 
in Fig. S3F. 

The BC pattern analysis can be extended to complex networks of neural cells. Let us 
consider the network of FHN model(l,2) 



which has been extensively used to describe neural dynamics. For the given parameter set 
given here the individual neural cells are excitable while nonoscillatory. For certain initial 
preparations the network of coupled cells can be self-organized to sustained oscillations. Fig. 
S4^ shows one of such network and Figs. S45 and 4 C present two different periodic orbits of 
the same network structure in Fig. S4^ for two sets of different initial conditions. Figs. SAD 
and S4E present the BC patterns corresponding to states Figs. S4B and S4C, respectively. 
The one circle (Fig. SAD) and two circle (Fig. S4£^ BC structures are identified by applying 
the CR method. 




dvj 
dt 



e{ui + /? - 7fi) 



i = l,2,...,Ar 



{SI) 
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Part 2 Oscillatory dynamics of genomic regulatory net- 
works of Eqs. (2) 

Firstly, the oscillatory negative feedback loops(NFLs) of GRNs for g ^ 5 are presented in 
Fig. S5. 

In order to define phase- advanced interactions of Eqs. (6) and (7), we should first specify 
the meaning of phase <f)i with Fourier decomposition. With (pi known, the phase of interaction 
from node j to node i (denoted by can be identified by (pj {(f) j+n) for positive (negative) 
interaction, and the interaction with < ~ 0i i$ f is called phase-advanced driving. 
For defining the dominant phase-advanced driving of a given node i, we should identify the 
interaction with the maximum amplitude among all the phase-advanced interactions of node 
i. 

The phase of single node i can be define as follows: 

V + Pi V «i + Pi 

ai = 7f j Axism{ — )dt, Pi = f J AXiCos[—)dt, 

_ 1 r 

Axi = Xi{t) - Xi, Xi = — Xi{t)dt, i = l,2,---,N. 

^ Jo 

Suppose Aj_,i{t) is the interaction from node j to node i which can be computed explicitly 
by linear approximation as 



AJ^^it) = ^Ax,, {S3) 



where ^ is a constant valued at time averages Xk {k = 1,2, ■ ■ ■ , N). Aj^i{t) is periodically 
oscillatory with period T and zero average. With Eqs. (S2) (S3) the phase (fyj^i can be define 
by (j)j {(j)j + n) for positive (negative) interaction, and the amplitude of Aj^i{t) is given by 
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The dominant phase-advanced driving of node i is the interaction with the largest ampHtude 
IIAj^ill among all the phase-advanced interactions of node i (i.e, all j s with < (pj^i — 4>i^ 
-) 

The crucial point of the approximations of Eqs. (S2)-(S4) is to neglect the contributions 
of all high-frequency harmonics of the interactions and thus the definitions work very well for 
the cases of single-peaked periodic motions. Some slight modification is needed if the motion 
is multiple-peaked, and these cases are not considered in the present article. It happens with 
extremely low probability that the sum of interactions of a given node obeys the driving 
condition Eq. (7) while individual interactions do not. In this case we simply define the 
interaction with phase nearest to the region as the dominant interaction. These cases occur 
always (with our observations) for the nodes near the branch ends. 



Part 3 Oscillatory dynamics of genomic regulatory net- 
works of type "OR" interactions 

Eqs. (2) consider "AND " type of joint regulations of activatory and repressive regulators, 
represented by the productive formula. In realistic regulatory networks there also exist joint 
regulations of type "OR ". Typical mathematical formulas of a network of this type (3,4) 
with size are as follows: 



no coupling , 

(1 - ^'i) j^iH^'^.hi positive couphng , 

k - ■ X 

-Xi ^. negative couphng . 

Both Eqs. (2) and (S5) arc approximations of more realistic as well as more complicated 
regulatory networks mixing "AND" and "OR" dynamics. 

Since in Eqs. (S5) each node does not oscillate individually and positive (negative) regu- 
lations tend to increase (reduce) monotonously the concentration of the protein interacted. 
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the phase relations analyzed in the article, are satisfied entirely. In particular, the phase 
and amplitude of couplings are given by Eqs. (6) and (S4), respectively, with fi replaced by 
^f=ifij- For mathematical simplicity we again take identical parameters in Eqs. (S5) for 
all nodes hi — h, kij — k, Ki — K, i — 1,2, N, except 7i=i,...,Ar-i — ^, — 1-0. 

In Figs. S6 we do the same as Fig. 5 with Eqs. (2) replaced by the "OR" dynamics of Eqs. 
(S5). It is really striking that though the coupling dynamics of Eqs. (S5) is different from 
that of Eqs. (2) considerably, BC patterns of Eqs. (S5) have essentially the same structures as 
Fig. 5, and all the features predicted in Fig. 2 are confirmed perfectly in Fig. S6. Moreover, 
the prediction of oscillation sources (one of the oscillatory NFLs in Fig. S5) is also fully 
verified by all the BC circles obtained from Eqs. (S5). We have also made statistics of the 
circuits of Eqs. (S5) with huge number of tests by varying coupling structures, parameter 
values and initial conditions (10® for 3-node, 10'^ for 4, 5-node, 10^ for 10, 20-node). We find 
that the predictions on BC structures in Fig. 2 and Fig. S5 are well confirmed by all tests 
of oscillatory circuits with no exception. 
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Fig. SI. Oscillatory behaviors of a complex ECN of Eqs. (1) with = 200, u = 3, 
Du — 1.0. All other parameters are the same as Fig. 1(A). (A) Interaction structure under 
consideration. (B)-(D) BC patterns of periodically oscillatory states of ECN (A) for three 
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different initial conditions. All red nodes have the same meanings as those in article. (B) 
BC pattern with a 27r circle. (C) BC pattern with two source circles. In order to suppress 
the oscillation one has to remove the two red nodes 6 and 38 simultaneously, destroying 
both source circles. (D) BC pattern with with a An (wave number k — 2) circle of circle 
size n — 21. We find two red nodes in the circle separated by about 27r distance. In order 
to suppress oscillation we have to remove both red nodes 126 and 178 simultaneously. (E) 
The BC pattern constructed for the new oscillatory state after removing 140 nodes from 
(D). Oscillation persists and the original BC circle remains unchanged after 70% nodes are 
removed. Similar results (not shown) can be obtained also for BC patterns of (B) and (C). 
Fig. S2. BC pattern manipulations. (A)(B) BC patterns after removing a red node 6 and 
the other red node 38 from Fig. S1(C), respectively. In both cases, the circle containing the 
removed node is destroyed and the remained circle serves as the unique oscillation source 
of the network. (C)(D) BC patterns with red nodes 126 and 178 removed from Fig. S1(D), 
respectively. Now the 47r BC circle of Fig. S1(D) shrinks to a 2tt circle with the remaining 
red node being the key circle node in the new BC pattern. 

Fig. S3. BC patterns of Eqs. (1) for different system sizes with u = 3, Du = 1.0. All other 
parameters are the same as Fig. 1(A). Interaction structures and initial conditions are chosen 
randomly. (A)-(C) BC patterns with A*" = 400 for different states. Numbers inside the circle 
indicate the corresponding indexes of corresponding nodes, and the numbers outside the 
circle, associated to circle nodes by arrows, denote the sizes of the branches rooted at the 
given nodes. All nodes without arrows have no downstream branch. Single 27r BC circle; 
single 47r circle and double-circle structure arc observed in (A)(B)(C), respectively. (D)-(F) 
The same as (A)-(C) with N — 2500. Single-, double- and triple-circle BC structures are 
identified in (D)-(F), respectively. 

Fig. S4. Oscillatory complex network of FHN cells (Eqs. (SI)) and some examples of BC 
patterns. The parameters are set as e = 0.2, 7 = 0.5, f3 = 1.0. All couphngs with coupling 
strength Du = 0.1 are randomly chosen between N nodes each having coupling degree u. 
(A) An example of FHN network with N — 100, u — 3. (B)(C) Two periodic orbits of 
network (A) for two different initial conditions. (D)(E) The BC patterns corresponding to 
the states (B) and (C), respectively. The red nodes have the same meanings as Figs. 3(A) 
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and 4(C). 

Fig. S5. All oscillatory g-node NFLs of Eqs. (2) for g < 5. The solid green (red dashed) line 
denotes positive (negative) regulation. (A) 2-node NFL. (B) 3-node NFL with three negative 
interactions. (C) 3-node NFL with a single negative interaction. (D) 4-node NFL with three 
negative interactions. (E) 4-node NFL with a single negative interaction. (F) 5- node NFL 
with a single negative interaction. (G) 5-node NFL with three negative interactions (Type I). 
(H) 5-node NFL with three negative interactions (Type II). (I) 5-node NFL with five negative 
interactions. All these oscillatory NFLs have been observed in various BC patterns. 
Fig. S6. Numerical results on oscillatory GRNs of Eqs. (S5) (type "OR"). Homogeneous 
parameter sets are used with hi — h, kij — k, Ki — K, i — 1,2, ... ,N, except 7i=i,...,Ar-i = 7, 
7Ar = 1.0. We use two sets of parameters in simulation. In Figs. (F), (G), (J) h = 3, 
k = 734.5525513, K = 3.5531352, 7 = 0.8991836; in Figs. (H), {!) h = A, k = 753.6924438, 
K — 0.6794546, 7 — 0.1966502. These parameters are randomly taken in the available ranges 
specified in Science report (3). In all systems unique periodically oscillatory solutions are 
asymptotically approached from arbitrary initial conditions. (A)-(E) Interaction structures 
of GRNs of different size N and different number / of interactions. (A) N = A, I = 9. (B) 
N = 6, I = 28. (C) N = 10, I = 56. (D) N = 20, I = 109. (E) N = AO, I = 222. (F)-(J) 
BC patterns corresponding to (A)-(E), respectively. Red and empty nodes have the same 
meanings as in Figs. 5 (F)-(J). 
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Fig. S3 



34 




35 




Fig. S5 
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Fig. S6 
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